% Woodford implemented as in Costain and Nakov, 2008
% permanent quality shocks
% ANTON NAKOV AND PETER KARADI (MARCH 2022)

clear all; clc; close all;

tic

global par_;

params_set;

% set initial GUESSES
switch model
    case 2  % Calvo 
            params0 = std_eA0;
    case 3  % Golosov-Lucas 
            params0 = [menucost0; std_eA0];
    case {4,5,7}  %Woodford, linear, uniform random mc
            params0 = log(alpha0);
    case 6  %asymmetric linear
            params0 = [log(alpha0); log(alphapos0)];      
end

%params_dyn0 = [tan(phiR0*(pi/2)); log(std_m0)];  
switch par_.switch_shock
    case 'm'
        params_dyn0 = [log(std_m0)];
    case 'A'
        params_dyn0 = [log(std_A0)];
end

% Pack parameters
par_.model = model;
par_.mymoment = mymoment;
par_.mfreq = mfreq;
par_.estim = estim;
par_.estim_type = estim_type;
par_.estim_dyn = estim_dyn;
par_.switch_sim = switch_sim;
par_.doss = doss;
par_.dodyn = dodyn; 
par_.show = false;
par_.comp = false;


opts_fsolve = optimset('Display','iter','TolFun',1e-6,'TolX',1e-6,'MaxIter',10);  
opts_fmin = optimset('Display','iter','TolX',sqrt(eps));  
opts_sa = optimset('Display','iter','TolFun',1e-6);  

par_.wbar0_min = wbar0_min;
par_.wbar0_max = wbar0_max;
par_.chi0_min  = chi0_min;
par_.chi0_max  = chi0_max;
par_.phiR      = phiR0;
par_.phiA      = phiA0;
par_.std_m     = std_m0;
par_.std_A     = std_A0;

if estim % find best fitting params in the steady state
     par_.dodyn = 0;
     par_.wbar   = wbar0;
     par_.chi    = chi0;
     par_.menucost = menucost0;
     par_.std_eA = std_eA0;
     par_.omega  = omega0;
     fprintf('\n\nEstimating model %0.1g \n\n',model) %#ok<*UNRCH>
     switch model
         case {2,3}
            [out,fval,flag] = fsolve('sdp_solve',params0,opts_fsolve);
         case {4,5,6,7}
                     %[out,fval,flag] = simulannealbnd(@sdp_solve,params0,-Inf,Inf,opts_fsolve);
                    %[out,fval,flag] = fminbnd('sdp_solve',1e-1,1e2,opts_fmin);
                    %[out,fval,flag] = simulannealbnd(@sdp_solve,params0,log(0.1),log(100),opts_sa);
                    switch estim_type
                        case 'quick' %needs good starting values, might stuck?
                            [out,fval,flag] = fsolve('sdp_solve',params0,opts_fsolve);
                        case 'robust' %need to be extended if used for model 6
                            %First systematically overview the range of relevant starting points
                            if exist('alpha_distance.mat','file')==0
                                alpha0_min = 0.01;
                                alpha0_max = 10;
                                no_alpha0   = 25;
                                logalpha0_inc  = (log(alpha0_max)-log(alpha0_min))/(no_alpha0-1);
                                logalpha0_vec = log(alpha0_min):logalpha0_inc:log(alpha0_max);
                                distance_vec = NaN(1,no_alpha0);
                                for ii=1:no_alpha0
                                    distance_vec(1,ii)=sdp_solve(logalpha0_vec(1,ii));
                                    save alpha_distance.mat logalpha0_vec distance_vec;
                                end
                                figure(1000); plot(exp(logalpha0_vec),distance_vec); xlabel('\alpha'); ylabel('Distance');                                
                            else
                                load alpha_distance.mat logalpha0_vec distance_vec;
                            end
                            distance_min_ind = find(distance_vec==min(distance_vec));
                            %Use the minimum as a starting value in a Newton algorithm
                            [out,fval,flag] = fsolve('sdp_solve',logalpha0_vec(1,distance_min_ind),opts_fsolve);
                    end
    end
     par_.show = true;
     sdp_solve_w(par_.wbar);
elseif estim_dyn  %estimate parameters related to the dynamics
     par_.dodyn = false;  
     par_.doss = true;
     par_.wbar   = wbar0;
     par_.chi    = chi0;
     par_.menucost = menucost0;
     par_.std_eA = std_eA0;
     par_.omega = omega0;
     sdp_solve(params0);   %solves for wbar, chi
     [out,fval,flag] = fsolve('fun_dynsim',params_dyn0,opts_fsolve);
     %par_.phiR = atan(out(1))/(pi/2);
     switch par_.switch_shock
         case 'm'
            par_.std_m = exp(out(1));
         case 'A'
             par_.std_A = exp(out(1));
     end
     par_.show = true;
     fun_dynsim(out);
elseif dodyn
     par_.dodyn = false;  
     par_.doss = true;
     par_.wbar   = wbar0;
     par_.chi    = chi0;
     par_.menucost = menucost0;
     par_.std_eA = std_eA0;
     par_.omega = omega0;
     sdp_solve(params0);   %solves for wbar, chi
     %[out,fval,flag] = fsolve('sdp_solve_chi',chi0,opts);
     par_.doss = false;
     par_.show = true;
     par_.dodyn = true;
     fun_dynsim(params_dyn0);
     %sdp_solve_w(par_.wbar);
elseif doss
     par_.dodyn = false;
     %[par_.chi,fval,flag] = fsolve('sdp_solve_chi',chi0,opts);
     par_.wbar   = wbar0;
     par_.chi    = chi0;
     par_.menucost = menucost0;
     par_.std_eA = std_eA0;
     par_.omega = omega0;
     sdp_solve(params0);   %solves for wbar, chi
     par_.show = true;
     sdp_solve_w(par_.wbar);
end
      
sound(sin(1:600)); toc